##############################
# STATE REACH AND DEVELOPMENT // PSRM // Replication
#
# Main Analysis File
#
##############################

# DATA #######################

# Load main data
data.ls <- list(educ = readRDS(file.path("data", "data_educ.rds")),
                infmort = readRDS(file.path("data", "data_infmort.rds")),
                nightlight = readRDS(file.path("data", "data_nightlightmain.rds")))


# SPECIFICATION ###############

# Data for panel analysis
main.data.names <- c("educ","infmort", "nightlight")


# Main endogenous / treatment variables
treat.vars <- c("natcap.ln" , "regcap.ln")
names(treat.vars) <-  c("Time to nat. capital (log)","Time to reg. capital (log)")

# Dependent variable
get_depvar <- function(data.name){
  if(data.name == "educ"){
    dep.var <- c(`Primary educ. (0/100)` = "educ.prim")
  } else if(data.name == "infmort"){
    dep.var <- c(`Infant mort. (0/100)` = "dead")
  } else if(data.name == "nightlight"){
    dep.var <- c(`Light/capita (log)` = "nightlight.pc.log")
  }
  return(dep.var)
}

# Dataset specific control variables
get_contr <- function(data.name){
  if(data.name == "educ"){
    contr <- c(`Female dummy` = "female" ,`Age` = "age.num"  , `Age sq.` = "I(age.num^2)")
  } else if(data.name == "infmort"){
    contr <-c(`Female dummy` = "female",`Birth order` = "birthorder.num" , "I(birthorder.num^2)",`Twin dummy` =  "twin" , 
              `Mother's age at birth` = "mother.b.age" , "I(mother.b.age^2)" )
  } else if(data.name == "ab"){
    contr <- c()
  } else if(data.name == "nightlight"){
    contr <- c()
  }
  return(contr)
}

# Fixed effects
fe.global <- "pts.id + cow.year + survey.id"

# Clustering of standard errors
cluster.global <- "pts.id + cow.year"


# Notes
latex.notes.childmort <- function(width = .95){
  paste0("\\parbox[t]{",width,"\\textwidth}{\\textit{Notes:} OLS linear models. 
         Control variables include an infant's mother's age and age squared, the birthorder and its square, 
         as well as a female and twin dummy. 
         Standard errors clustered on the point and country-year levels.
         Significance codes: $^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01}")}
latex.notes.educ <-  function(width = .95){
  paste0("\\parbox[t]{",width,"\\textwidth}{\\textit{Notes:} OLS linear models. 
         Control variables include respondents' age and age squared, as well as a female dummy. 
         Standard errors clustered on the point and country-year levels.
         Significance codes: $^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01}")}
latex.notes.nl <-  function(width = .95){
  paste0("\\parbox[t]{",width,"\\textwidth}{\\textit{Notes:} OLS linear models. 
         Standard errors  clustered on the point and country-year levels.
         Significance codes: $^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01}")}
latex.notes.comb <-  function(width = .95){
  paste0("\\parbox[t]{",width,"\\textwidth}{\\textit{Notes:} OLS linear models. 
  Individuals are the units of the education and infant mortality analyses, Voronoi cells those of the nightlight analyses. 
         Control variables for models with primary education as the dependent variable consist of
         responents' age and age squared, as well as a female dummy. Where infant mortality is the dependent variable,
         models include an infant's mother's age at birth and its square, the birthorder and its square, 
         as well as a female and twin dummy. 
         Standard errors clustered on the point and country-year levels.
         Significance codes: $^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01}")}


# Coefficient list
all.coef.list <- list()



##################################
# DESCRIPTIVES ###################
##################################

source("code/descriptives.R")



##################################
# TABLE 1: MAIN SPECIFICATION ####
##################################

stub <- "tab1_main"

# ... estimate
model.ls <- lapply(main.data.names, function(dn){
  form <- make_form(dv = get_depvar(dn), expl = c(treat.vars,get_contr(dn)), 
                    fe = fe.global, 
                    iv = "0",
                    se = cluster.global )
  m <- felm(form,  data = data.ls[[dn]], keepX = F, keepCX = F)
  m
})

#... sum of betas
sum.betas <- lapply(model.ls, function(m){
  dm <- deltaMethod(m, paste(treat.vars, collapse = " + "), vcov. = m$clustervcv)
  wt <- lfe::waldtest(m, as.formula(paste("~", paste(treat.vars, collapse = " + "))), type = "cluster")
  dm$pvalue <- wt["p"]
  dm
})

# ... prepare
add.lines <- list(latex.addline_nomc("$\\beta_1 + \\beta_2$:", 
                                     unlist(lapply(sum.betas, function(x){
                                       paste0(round(x[["Estimate"]], 3), 
                                              ifelse(x[["pvalue"]] < .01, "^{***}",ifelse(x[["pvalue"]] < .05, "^{**}",ifelse(x[["pvalue"]] < .1, "^{*}","" ) ) ))
                                     }))),
                  latex.addline_nomc("", unlist(lapply(sum.betas, function(x) paste0("(",round(x[["SE"]], 3), ")")))),
                  latex.addline("Point FE:", rep("yes", length(model.ls))),
                  latex.addline("Country-year FE:", rep("yes", length(model.ls))),
                  latex.addline("Survey FE:", rep(c("yes","yes","--"), 1)),
                  latex.addline("Controls:", rep(c("yes","yes","--"), 1)),
                  latex.mean.dv(model.ls))
add.lines[[2]][[length(add.lines[[2]])]] <- paste(add.lines[[2]][[length(add.lines[[2]])]], "\\\\ \\hline \\\\[-1.8ex] ")

# ... save main
fileConn <- file(file.path(tab.path, paste0(stub,".tex")))
writeLines(stargazer(model.ls,
                     title="Changes in time to national/regional capital and local development",
                     keep = treat.vars,
                     multicolumn = F,
                     column.labels = names(unlist(lapply(main.data.names, get_depvar))),
                     dep.var.caption = "",dep.var.labels = rep("", length(model.ls)),
                     covariate.labels = names(treat.vars),
                     font.size = "scriptsize",
                     notes.align = "c", label=stub, align =T,
                     add.lines = add.lines,digits = 3, intercept.top = T,intercept.bottom = F,
                     omit.stat = c("rsq","res.dev","ser"),
                     notes = latex.notes.comb(.85), notes.label = "", notes.append = F), 
           fileConn)
close(fileConn)

# Save coefs
all.coef.list[[stub]] <- extract_coef(model.ls)

# TABLE A6: NATIONAL CAPITAL DROP AND KEEP OUTLIER CASES ###
stub <- "taba6_dropkeepcases"

# ... estimate
model.ls <- unlist(lapply(main.data.names, function(dn){
  if(dn %in% c("educ", "infmort")){
    cases <- c(437, 475, 510)
  } else {
    cases <- c(626, 531)
  }
  
  lapply(c(0, 1), function(c){
    print(c)
    form <- make_form(dv = get_depvar(dn), expl = c(treat.vars,get_contr(dn)), 
                      fe = fe.global, 
                      iv = "0",
                      se = cluster.global )

    if(dn == "nightlight"){

      ## make stable cowcodes by correcting Sudan and Eritrea
      data <- data.ls[[dn]]
      data$cowcode.stable <- data$cowcode
      data$cowcode.stable[data$pts.id %in% unique(data$pts.id[data$cowcode == 626])] <- 626
      data$cowcode.stable[data$pts.id %in% unique(data$pts.id[data$cowcode == 531])] <- 531

      ## Estimate models
      if(c == 0){
        m <- felm(form,  data = data[!data$cowcode.stable %in% cases,], keepX = F, keepCX = F)
      } else {
        m <- felm(form,  data = data[data$cowcode.stable %in% cases,], keepX = F, keepCX = F)
      }
    } else {
      if(c == 0){
        m <- felm(form,  data = data.ls[[dn]][!data.ls[[dn]]$cowcode %in% cases,], keepX = F, keepCX = F)
      } else {
        m <- felm(form,  data = data.ls[[dn]][data.ls[[dn]]$cowcode %in% cases,], keepX = F, keepCX = F)
      }
    }
    
  })
}), recursive = F)

#... sum of betas
sum.betas <- lapply(model.ls, function(m){
  dm <- deltaMethod(m, paste(treat.vars, collapse = " + "), vcov. = m$clustervcv)
  wt <- lfe::waldtest(m, as.formula(paste("~", paste(treat.vars, collapse = " + "))), type = "cluster")
  dm$pvalue <- wt["p"]
  dm
})

# ... prepare
add.lines <- list(latex.addline_nomc("$\\beta_1 + \\beta_2$:", 
                                     unlist(lapply(sum.betas, function(x){
                                       paste0(round(x[["Estimate"]], 3), 
                                              ifelse(x[["pvalue"]] < .01, "^{***}",ifelse(x[["pvalue"]] < .05, "^{**}",ifelse(x[["pvalue"]] < .1, "^{*}","" ) ) ))
                                     }))),
                  latex.addline_nomc("", unlist(lapply(sum.betas, function(x) paste0("(",round(x[["SE"]], 3), ")")))),
                  latex.addline("Cap.-re-loc. cases:", rep(c("drop","only"), length(main.data.names))),
                  latex.addline("Point FE:", rep("yes", length(model.ls))),
                  latex.addline("Country-year FE:", rep("yes", length(model.ls))),
                  latex.addline("Survey FE:", rep(c("yes","yes","--"), each = 2)),
                  latex.addline("Controls:", rep(c("yes","yes","--"), each = 2)),
                  latex.mean.dv(model.ls))
add.lines[[2]][[length(add.lines[[2]])]] <- paste(add.lines[[2]][[length(add.lines[[2]])]], "\\\\ \\hline \\\\[-1.8ex] ")

# ... save main
fileConn <- file(file.path(tab.path, paste0(stub,".tex")))
writeLines(stargazer(model.ls,
                     title="Changes in time to national/regional capital and local development: Cases where national capitals changed",
                     keep = treat.vars,
                     multicolumn = F,
                     column.labels =collab_w_ruler(column.labels = names(unlist(lapply(main.data.names, get_depvar))),
                                                   column.separate = c(2,2,2)),
                     column.separate = c(2,2,2), column.sep.width ="-10pt",
                     dep.var.caption = "",dep.var.labels = rep("", length(model.ls)),
                     covariate.labels = names(treat.vars),
                     font.size = "scriptsize",
                     notes.align = "c", label=stub, align =T,
                     add.lines = add.lines,digits = 3, intercept.top = T,intercept.bottom = F,
                     omit.stat = c("rsq","res.dev","ser"),
                     notes = latex.notes.comb(), notes.label = "", notes.append = F), 
           fileConn)
close(fileConn)

# Save coefs
all.coef.list[[stub]] <- extract_coef(model.ls)



# TABLE A7: CONTROL FOR MIGRATION ##########
stub <- "taba7_migration"

# Load individual recode data with migration variable
ir.mig.df <- readRDS(file.path("data", "data_educ_mig.rds"))

# Estimate
model.ls <- list(felm(make_form(dv = get_depvar("educ"), 
                                expl = c(treat.vars, "is.migrant", paste0("I(is.migrant*", treat.vars, ")"), 
                                         get_contr("educ")), 
                                fe = fe.global, 
                                iv = "0",
                                se = cluster.global ),  
                      data = ir.mig.df, keepX = F, keepCX = F),
                 felm(make_form(dv = get_depvar("educ"), 
                                expl = c(treat.vars, "is.migrant", paste0("I(is.migrant*", treat.vars, ")"), 
                                         paste0("is.migrant*", get_contr("educ"))), 
                                fe = fe.global, 
                                iv = "0",
                                se = cluster.global ),  
                      data = ir.mig.df, keepX = F, keepCX = F),
                 felm(make_form(dv = get_depvar("infmort"), 
                                expl = c(treat.vars, "is.migrant", paste0("I(is.migrant*", treat.vars, ")"), 
                                         get_contr("infmort")), 
                                fe = fe.global, 
                                iv = "0",
                                se = cluster.global ),  
                      data = data.ls[["infmort"]], keepX = F, keepCX = F),
                 felm(make_form(dv = get_depvar("infmort"), 
                                expl = c(treat.vars, "is.migrant", paste0("I(is.migrant*", treat.vars, ")"), 
                                         paste0("is.migrant*",get_contr("educ"))), 
                                fe = fe.global, 
                                iv = "0",
                                se = cluster.global ),  
                      data = data.ls[["infmort"]], keepX = F, keepCX = F))


#... sum of betas
sum.betas <- lapply(model.ls, function(m){
  dm <- deltaMethod(m, paste(treat.vars, collapse = " + "), vcov. = m$clustervcv)
  wt <- lfe::waldtest(m, as.formula(paste("~", paste(treat.vars, collapse = " + "))), type = "cluster")
  dm$pvalue <- wt["p"]
  dm
})

# ... prepare
add.lines <- list(latex.addline_nomc("Non-migrants: $\\beta_1 + \\beta_2$:", 
                                     unlist(lapply(sum.betas, function(x){
                                       paste0(round(x[["Estimate"]], 3), 
                                              ifelse(x[["pvalue"]] < .01, "^{***}",ifelse(x[["pvalue"]] < .05, "^{**}",ifelse(x[["pvalue"]] < .1, "^{*}","" ) ) ))
                                     }))),
                  latex.addline_nomc("", unlist(lapply(sum.betas, function(x) paste0("(",round(x[["SE"]], 3), ")")))),
                  latex.addline("Migrant $\\times$ controls", c("no", "yes", "no", "yes")),
                  latex.addline("Point FE:", rep("yes", length(model.ls))),
                  latex.addline("Country-year FE:", rep("yes", length(model.ls))),
                  latex.addline("Survey FE:", rep(c("yes","yes","--"), each = 2)),
                  latex.addline("Controls:", rep(c("yes","yes","--"), each = 2)),
                  latex.mean.dv(model.ls))
add.lines[[2]][[length(add.lines[[2]])]] <- paste(add.lines[[2]][[length(add.lines[[2]])]], "\\\\ \\hline \\\\[-1.8ex] ")

# ... save main
fileConn <- file(file.path(tab.path, paste0(stub,".tex")))
writeLines(stargazer(model.ls,
                     title="Changes in time to national/regional capital and local development: Migrants and non-migrants",
                     keep = 1:5,
                     multicolumn = F,# se = se,
                     column.labels =collab_w_ruler(column.labels = names(unlist(lapply(c("educ","infmort"), get_depvar))), 
                                                   column.separate = rep(2, 2), 
                                                   trim = 10, add.below = NULL),
                     column.separate = rep(2, 2), 
                     column.sep.width ="0pt",
                     dep.var.caption = "",dep.var.labels = rep("", length(model.ls)),
                     covariate.labels = c(names(treat.vars), "Migrant", paste0("Migrant$\\times$", names(treat.vars))),
                     font.size = "scriptsize",
                     notes.align = "c", label=stub, align =T,
                     add.lines = add.lines,digits = 3, intercept.top = T,intercept.bottom = F,
                     omit.stat = c("rsq","res.dev","ser"),
                     notes = latex.notes.comb(), notes.label = "", notes.append = F), 
           fileConn)
close(fileConn)

# Save coefs
all.coef.list[[stub]] <- extract_coef(model.ls)


# TABLES A8 AND A9: LEAD TREATMENT INDICATORS ######
source("code/analysis_leads.R")




# TABLE A10: CONTROL FOR MARKET ACCESS ######
stub <- "taba10_maaccess"

# Prepare
these.add.vars <- paste0("log(.001+",c("ma.int.donal - ma.nat.donal" ,  
                                       "ma.int.eaton - ma.nat.eaton", 
                                       "ma.nat.donal", "ma.nat.eaton"),")")

# Estimate
model.ls <- lapply(main.data.names, function(dn){
  form <- make_form(dv = get_depvar(dn), 
                    expl = c(treat.vars, these.add.vars, 
                             get_contr(dn)), 
                    fe = fe.global, 
                    iv = "0",
                    se = cluster.global )
  felm(form,  data = data.ls[[dn]], keepX = F, keepCX = F)
})


#... sum of betas
sum.betas <- lapply(model.ls, function(m){
  dm <- deltaMethod(m, paste(treat.vars, collapse = " + "), vcov. = m$clustervcv)
  wt <- lfe::waldtest(m, as.formula(paste("~", paste(treat.vars, collapse = " + "))), type = "cluster")
  dm$pvalue <- wt["p"]
  dm
})

# ... prepare
add.lines <- list(latex.addline_nomc("$\\beta_1 + \\beta_2$:", 
                                     unlist(lapply(sum.betas, function(x){
                                       paste0(round(x[["Estimate"]], 3), 
                                              ifelse(x[["pvalue"]] < .01, "^{***}",ifelse(x[["pvalue"]] < .05, "^{**}",ifelse(x[["pvalue"]] < .1, "^{*}","" ) ) ))
                                     }))),
                  latex.addline_nomc("", unlist(lapply(sum.betas, function(x) paste0("(",round(x[["SE"]], 3), ")")))),
                  latex.addline("Point FE:", rep("yes", length(model.ls))),
                  latex.addline("Country-year FE:", rep("yes", length(model.ls))),
                  latex.addline("Survey FE:", rep(c("yes","yes","--"), 1)),
                  latex.addline("Controls:", rep(c("yes","yes","--"), 1)),
                  latex.mean.dv(model.ls))
add.lines[[2]][[length(add.lines[[2]])]] <- paste(add.lines[[2]][[length(add.lines[[2]])]], "\\\\ \\hline \\\\[-1.8ex] ")

# ... save main
fileConn <- file(file.path(tab.path, paste0(stub,".tex")))
writeLines(stargazer(model.ls,
                     title="Changes in time to national/regional capital and local development: Controlling for market access",
                     keep = c(treat.vars,"ma.int.donal", "ma.nat.donal" ,  "ma.int.eaton", "ma.nat.eaton"),
                     multicolumn = F,
                     column.labels = names(unlist(lapply(main.data.names, get_depvar))),
                     column.separate = rep(length(these.add.vars), length(main.data.names)), 
                     column.sep.width ="5pt",
                     dep.var.caption = "",dep.var.labels = rep("", length(model.ls)),
                     covariate.labels = c(names(treat.vars), 
                                          "MA, internat. (log; $\\theta=3.8$)","MA, nat. (log; $\\theta=3.8$)",
                                          "MA, internat. (log; $\\theta=8.28$)","MA, nat. (log; $\\theta=8.28$)"),
                     font.size = "scriptsize",
                     notes.align = "c", label=stub, align =T,
                     add.lines = add.lines,digits = 3, intercept.top = T,intercept.bottom = F,
                     omit.stat = c("rsq","res.dev","ser"),
                     notes = latex.notes.comb(.85), notes.label = "", notes.append = F), 
           fileConn)
close(fileConn)

# Save coefs
all.coef.list[[stub]] <- extract_coef(model.ls)



# TABLE A11: CONTROL TIME-VARYING ETHNIC POLITICS #####
stub <- "taba11_ethnicpolitics"

# Prepare
these.add.vars <- c("eprincl","eprexcl","ecwar","I(cumsum.ecwar>0)",
                    "timesince.ecwar","I(timesince.ecwar^2)")

# Estimate
model.ls <- lapply(main.data.names, function(dn){
  form <- make_form(dv = get_depvar(dn), 
                    expl = c(treat.vars, these.add.vars,get_contr(dn)), 
                    fe = fe.global, 
                    iv = "0",
                    se = cluster.global )
  felm(form,  data = data.ls[[dn]], keepX = F, keepCX = F)
})

#... sum of betas
sum.betas <- lapply(model.ls, function(m){
  dm <- deltaMethod(m, paste(treat.vars, collapse = " + "), vcov. = m$clustervcv)
  wt <- lfe::waldtest(m, as.formula(paste("~", paste(treat.vars, collapse = " + "))), type = "cluster")
  dm$pvalue <- wt["p"]
  dm
})

# ... prepare
add.lines <- list(latex.addline_nomc("$\\beta_1 + \\beta_2$:", 
                                     unlist(lapply(sum.betas, function(x){
                                       paste0(round(x[["Estimate"]], 3), 
                                              ifelse(x[["pvalue"]] < .01, "^{***}",ifelse(x[["pvalue"]] < .05, "^{**}",ifelse(x[["pvalue"]] < .1, "^{*}","" ) ) ))
                                     }))),
                  latex.addline_nomc("", unlist(lapply(sum.betas, function(x) paste0("(",round(x[["SE"]], 3), ")")))),
                  latex.addline("Point FE:", rep("yes", length(model.ls))),
                  latex.addline("Country-year FE:", rep("yes", length(model.ls))),
                  latex.addline("Survey FE:", rep(c("yes","yes","--"), 1)),
                  latex.addline("Controls:", rep(c("yes","yes","--"), 1)),
                  latex.mean.dv(model.ls))
add.lines[[2]][[length(add.lines[[2]])]] <- paste(add.lines[[2]][[length(add.lines[[2]])]], "\\\\ \\hline \\\\[-1.8ex] ")

# ... save main
fileConn <- file(file.path(tab.path, paste0(stub,".tex")))
writeLines(stargazer(model.ls,
                     title="Changes in time to national/regional capital and local development: Controlling for ethnic representation and civil war",
                     keep = c(treat.vars, "eprincl","eprexcl","ecwar","cumsum.ecwar","timesince.ecwar"),
                     multicolumn = F,
                     column.labels = names(unlist(lapply(main.data.names, get_depvar))),
                     column.separate = rep(length(these.add.vars), length(main.data.names)), 
                     column.sep.width ="5pt",
                     dep.var.caption = "",dep.var.labels = rep("", length(model.ls)),
                     covariate.labels = c(names(treat.vars),"Ethnic inclusion (0/1)",  "Ethnic exclusion (0/1)",
                                          "Ethnic civil war (0/1)","Eth. war since indep. (0/1)", 
                                          "Time since eth. war", "Time since eth. war$^2$"),
                     font.size = "scriptsize",
                     notes.align = "c", label=stub, align =T,
                     add.lines = add.lines,digits = 3, intercept.top = T,intercept.bottom = F,
                     omit.stat = c("rsq","res.dev","ser"),
                     notes = latex.notes.comb(), notes.label = "", notes.append = F), 
           fileConn)
close(fileConn)

# Save coefs
all.coef.list[[stub]] <- extract_coef(model.ls)



# TABLE A12: CONTROLLING FOR CAPITALS ########

stub <- "taba12_capitalcontr"

# Prepare
add.capital.vars <- list(c("is.natcap", "is.regcap"),
                         c("I(natcap.ln < log(2))", "I(regcap.ln < log(2))"))

# Estimate
model.ls <- unlist(lapply(main.data.names, function(dn){
  lapply(add.capital.vars, function(av){
    form <- make_form(dv = get_depvar(dn), expl = c(treat.vars, av,get_contr(dn)), 
                      fe = fe.global, 
                      iv = "0",
                      se = cluster.global )
    felm(form,  data = data.ls[[dn]], keepX = F, keepCX = F)
  })
}), recursive = F)


#... sum of betas
sum.betas <- lapply(model.ls, function(m){
  dm <- deltaMethod(m, paste(treat.vars, collapse = " + "), vcov. = m$clustervcv)
  wt <- lfe::waldtest(m, as.formula(paste("~", paste(treat.vars, collapse = " + "))), type = "cluster")
  dm$pvalue <- wt["p"]
  dm
})

# ... prepare
add.lines <- list(latex.addline_nomc("$\\beta_1 + \\beta_2$:", 
                                     unlist(lapply(sum.betas, function(x){
                                       paste0(round(x[["Estimate"]], 3), 
                                              ifelse(x[["pvalue"]] < .01, "^{***}",ifelse(x[["pvalue"]] < .05, "^{**}",ifelse(x[["pvalue"]] < .1, "^{*}","" ) ) ))
                                     }))),
                  latex.addline_nomc("", unlist(lapply(sum.betas, function(x) paste0("(",round(x[["SE"]], 3), ")")))),
                  latex.addline("Point FE:", rep("yes", length(model.ls))),
                  latex.addline("Country-year FE:", rep("yes", length(model.ls))),
                  latex.addline("Survey FE:", rep(c("yes","yes","--"), each = 2)),
                  latex.addline("Controls:", rep(c("yes","yes","--"), each = 2)),
                  latex.mean.dv(model.ls))
add.lines[[2]][[length(add.lines[[2]])]] <- paste(add.lines[[2]][[length(add.lines[[2]])]], "\\\\ \\hline \\\\[-1.8ex] ")

# ... save main
fileConn <- file(file.path(tab.path, paste0(stub,".tex")))
writeLines(stargazer(model.ls,
                     title="Changes in time to national/regional capital and local development: Controlling for capitals",
                     keep = 1:6,
                     multicolumn = F,# se = se,
                     column.labels = collab_w_ruler(column.labels = names(unlist(lapply(main.data.names, get_depvar))), 
                                                    column.separate = rep(length(add.capital.vars), length(main.data.names)), 
                                                    trim = 14, add.below = NULL),
                     column.separate = rep(length(add.capital.vars), length(main.data.names)), 
                     column.sep.width ="-13pt",
                     dep.var.caption = "",dep.var.labels = rep("", length(model.ls)),
                     covariate.labels = c(names(treat.vars), "National capital (0/1)", "Regional capital (0/1)",
                                          "Time to nat. cap. $<$ 1hr","Time to reg. cap $<$ 1hr"),
                     font.size = "scriptsize",
                     notes.align = "c", label=stub, align =T,
                     add.lines = add.lines,digits = 3, intercept.top = T,intercept.bottom = F,
                     omit.stat = c("rsq","res.dev","ser"),
                     notes = latex.notes.comb(), notes.label = "", notes.append = F), 
           fileConn)
close(fileConn)

# Save coefs
all.coef.list[[stub]] <- extract_coef(model.ls)




# TABLE A13: ALTERNATIVE EDUCATION OUTCOMES #############
stub <- "taba13_alteducvars"


# Prepare
educ.outcomes <- c(`Any educ. (0/100)` = "I(100*(educ_years > 0))",`Educ. years (linear)` = "educ_years",
                   `Educ. years (log)` = "log(1+educ_years)", `Sec. educ. (0/100)` = "I((educ.num > 1)*100)")

# Code implausibly large values // missing as missing
data.ls[["educ"]]$educ_years <- as.numeric(data.ls[["educ"]]$educ_years)
data.ls[["educ"]]$educ_years[data.ls[["educ"]]$educ_years > 30] <- NA

# Estimate
model.ls <- c(lapply(educ.outcomes[1:3], function(dv){
  print(dv)
  felm(make_form(dv = dv, expl = c(treat.vars,get_contr("educ")), 
                 fe = fe.global, 
                 iv = "0",
                 se = cluster.global ),
       data = data.ls[["educ"]])
}),lapply(educ.outcomes[4:4], function(dv){
  print(dv)
  felm(make_form(dv = dv, expl = c(treat.vars,get_contr("educ")), 
                 fe = fe.global, 
                 iv = "0",
                 se = cluster.global ),
       data = data.ls[["educ"]][data.ls[["educ"]]$age.num >= 18,])
}))


#... sum of betas
sum.betas <- lapply(model.ls, function(m){
  dm <- deltaMethod(m, paste(treat.vars, collapse = " + "), vcov. = m$clustervcv)
  wt <- lfe::waldtest(m, as.formula(paste("~", paste(treat.vars, collapse = " + "))), type = "cluster")
  dm$pvalue <- wt["p"]
  dm
})

# ... prepare
add.lines <- list(latex.addline_nomc("$\\beta_1 + \\beta_2$:", 
                                     unlist(lapply(sum.betas, function(x){
                                       paste0(round(x[["Estimate"]], 3), 
                                              ifelse(x[["pvalue"]] < .01, "^{***}",ifelse(x[["pvalue"]] < .05, "^{**}",ifelse(x[["pvalue"]] < .1, "^{*}","" ) ) ))
                                     }))),
                  latex.addline_nomc("", unlist(lapply(sum.betas, function(x) paste0("(",round(x[["SE"]], 3), ")")))),
                  latex.addline("Point FE:", rep("yes", length(model.ls))),
                  latex.addline("Country-year FE:", rep("yes", length(model.ls))),
                  latex.addline("Survey FE:", rep("yes", length(model.ls))),
                  latex.addline("Controls:", rep(c("yes"), each = 4)),
                  latex.mean.dv(model.ls))
add.lines[[2]][[length(add.lines[[2]])]] <- paste(add.lines[[2]][[length(add.lines[[2]])]], "\\\\ \\hline \\\\[-1.8ex] ")

# ... save main
fileConn <- file(file.path(tab.path, paste0(stub,".tex")))
writeLines(stargazer(model.ls,
                     title="Changes in time to national/regional capital and various education outcomes",
                     keep = c(treat.vars),
                     multicolumn = F,# se = se,
                     column.labels = names(educ.outcomes),
                     column.separate = rep(1, length(model.ls)), 
                     column.sep.width ="-1pt",
                     dep.var.caption = "",dep.var.labels = rep("", length(model.ls)),
                     covariate.labels = c(names(treat.vars)),
                     font.size = "scriptsize",
                     notes.align = "c", label=stub, align =T,
                     add.lines = add.lines,digits = 3, intercept.top = T,intercept.bottom = F,
                     omit.stat = c("rsq","res.dev","ser"),
                     notes = latex.notes.educ(), notes.label = "", notes.append = F), 
           fileConn)
close(fileConn)

# Save coefs
all.coef.list[[stub]] <- extract_coef(model.ls)



# TABLE A14: ALTERNATIVE CHILD HEALTHCARE OUTCOMES ######
stub <- "taba14_healthqual"

# Outcomes
hc.qual.out <- c("any_prenatal", "public.birth", "any_assist",  "I(assist_traditional*100)")
names(hc.qual.out) <- c( "Prof. prenatal care", "Birth in public inst.", 
                         "Prof. birth assist.",  "Trad. birth assist.")
# Estimate
model.ls <- lapply(hc.qual.out, function(dv){
  print(dv)
  felm(make_form(dv = dv, expl = c(treat.vars,get_contr("infmort")), 
                 fe = fe.global, 
                 iv = "0",
                 se = cluster.global ),
       data = data.ls[["infmort"]])
})


#... sum of betas
sum.betas <- lapply(model.ls, function(m){
  dm <- deltaMethod(m, paste(treat.vars, collapse = " + "), vcov. = m$clustervcv)
  wt <- lfe::waldtest(m, as.formula(paste("~", paste(treat.vars, collapse = " + "))), type = "cluster")
  dm$pvalue <- wt["p"]
  dm
})

# ... prepare
add.lines <- list(latex.addline_nomc("$\\beta_1 + \\beta_2$:", 
                                     unlist(lapply(sum.betas, function(x){
                                       paste0(round(x[["Estimate"]], 3), 
                                              ifelse(x[["pvalue"]] < .01, "^{***}",ifelse(x[["pvalue"]] < .05, "^{**}",ifelse(x[["pvalue"]] < .1, "^{*}","" ) ) ))
                                     }))),
                  latex.addline_nomc("", unlist(lapply(sum.betas, function(x) paste0("(",round(x[["SE"]], 3), ")")))),
                  latex.addline("Point FE:", rep("yes", length(model.ls))),
                  latex.addline("Country-year FE:", rep("yes", length(model.ls))),
                  latex.addline("Survey FE:", rep("yes", length(model.ls))),
                  latex.addline("Controls:", rep(c("yes"), each = 4)),
                  latex.mean.dv(model.ls))
add.lines[[2]][[length(add.lines[[2]])]] <- paste(add.lines[[2]][[length(add.lines[[2]])]], "\\\\ \\hline \\\\[-1.8ex] ")

# ... save main
fileConn <- file(file.path(tab.path, paste0(stub,".tex")))
writeLines(stargazer(model.ls,
                     title="Changes in time to national/regional capital and quality of prenatal care and birth assistance (in percent)",
                     keep = c(treat.vars),
                     multicolumn = F,# se = se,
                     column.labels = names(hc.qual.out),
                     column.separate = rep(1, length(model.ls)), 
                     column.sep.width ="-1pt",
                     dep.var.caption = "",dep.var.labels = rep("", length(model.ls)),
                     covariate.labels = c(names(treat.vars)),
                     font.size = "scriptsize",
                     notes.align = "c", label=stub, align =T,
                     add.lines = add.lines,digits = 3, intercept.top = T,intercept.bottom = F,
                     omit.stat = c("rsq","res.dev","ser"),
                     notes = latex.notes.childmort(), notes.label = "", notes.append = F), 
           fileConn)
close(fileConn)

# Save coefs
all.coef.list[[stub]] <- extract_coef(model.ls)


# TABLE A15: ALTERNATIVE NIGHTLIGHT OUTCOMES ######
stub <- "taba15_nloutcomes"

# Prepare
outcomes <- c( "nightlight.pc.log", "nightlight.log", "nightlight.any")
names(outcomes) <- c( "Light/capita (log)", "Light (log)","Any Light (0/100)")

# Estimate
model.ls <- lapply(outcomes, function(dv){
  print(dv)
  felm(make_form(dv = dv, expl = c(treat.vars,get_contr("nightlight")), 
                 fe = fe.global, 
                 iv = "0",
                 se = cluster.global ),
       data = data.ls[["nightlight"]])
})


#... sum of betas
sum.betas <- lapply(model.ls, function(m){
  dm <- deltaMethod(m, paste(treat.vars, collapse = " + "), vcov. = m$clustervcv)
  wt <- lfe::waldtest(m, as.formula(paste("~", paste(treat.vars, collapse = " + "))), type = "cluster")
  dm$pvalue <- wt["p"]
  dm
})

# ... prepare
add.lines <- list(latex.addline_nomc("$\\beta_1 + \\beta_2$:", 
                                     unlist(lapply(sum.betas, function(x){
                                       paste0(round(x[["Estimate"]], 3), 
                                              ifelse(x[["pvalue"]] < .01, "^{***}",ifelse(x[["pvalue"]] < .05, "^{**}",ifelse(x[["pvalue"]] < .1, "^{*}","" ) ) ))
                                     }))),
                  latex.addline_nomc("", unlist(lapply(sum.betas, function(x) paste0("(",round(x[["SE"]], 3), ")")))),
                  latex.addline("Unit FE:", rep("yes", length(model.ls))),
                  latex.addline("Country-year FE:", rep("yes", length(model.ls))),
                  latex.addline("Controls:", rep(c("--"), each = 3)),
                  latex.mean.dv(model.ls))
add.lines[[2]][[length(add.lines[[2]])]] <- paste(add.lines[[2]][[length(add.lines[[2]])]], "\\\\ \\hline \\\\[-1.8ex] ")

# ... save main
fileConn <- file(file.path(tab.path, paste0(stub,".tex")))
writeLines(stargazer(model.ls,
                     title="Changes in time to national/regional capital and various nightlight measures",
                     keep = c(treat.vars),
                     multicolumn = F,# se = se,
                     column.labels = names(outcomes),
                     column.separate = rep(1, length(model.ls)), 
                     column.sep.width ="5pt",
                     dep.var.caption = "",dep.var.labels = rep("", length(model.ls)),
                     covariate.labels = c(names(treat.vars)),
                     font.size = "scriptsize",
                     notes.align = "c", label=stub, align =T,
                     add.lines = add.lines,digits = 3, intercept.top = T,intercept.bottom = F,
                     omit.stat = c("rsq","res.dev","ser"),
                     notes = latex.notes.nl(.85), notes.label = "", notes.append = F), 
           fileConn)
close(fileConn)

# Save coefs
all.coef.list[[stub]] <- extract_coef(model.ls)



# TABLE A16: COUNTRY / COHORT WEIGHTS #######

stub <- "taba16_weights"

# Estimate
model.ls <- unlist(lapply(main.data.names, function(dn){
  if(dn == "nightlight"){
    weight.vars <- c("population")
  } else {
    weight.vars <- c("cowcode", "cow.year")
  }
  lapply(weight.vars, function(w){
    print(w)
    form <- make_form(dv = get_depvar(dn), expl = c(treat.vars,get_contr(dn)), 
                      fe = fe.global, 
                      iv = "0",
                      se = cluster.global )
    this.data <- data.ls[[dn]]
    this.data <- na.omit(this.data[,unique(c(w, all.vars(form)))])
    if(w %in% c("cowcode", "cow.year")){
      weights <- gen_weights(this.data[,w])
    } else {  
      weights <- this.data[,w]
    }
    m <- felm(form, data = this.data, weights = weights, keepX = F, keepCX = F)
    m
  })
}), recursive = F)


#... sum of betas
sum.betas <- lapply(model.ls, function(m){
  dm <- deltaMethod(m, paste(treat.vars, collapse = " + "), vcov. = m$clustervcv)
  wt <- lfe::waldtest(m, as.formula(paste("~", paste(treat.vars, collapse = " + "))), type = "cluster")
  dm$pvalue <- wt["p"]
  dm
})

# ... prepare
add.lines <- list(latex.addline_nomc("$\\beta_1 + \\beta_2$:", 
                                     unlist(lapply(sum.betas, function(x){
                                       paste0(round(x[["Estimate"]], 3), 
                                              ifelse(x[["pvalue"]] < .01, "^{***}",ifelse(x[["pvalue"]] < .05, "^{**}",ifelse(x[["pvalue"]] < .1, "^{*}","" ) ) ))
                                     }))),
                  latex.addline_nomc("", unlist(lapply(sum.betas, function(x) paste0("(",round(x[["SE"]], 3), ")")))),
                  latex.addline("Weights", c(rep(c("country", "country"), 2), "population")),
                  latex.addline("", rep(c("", "cohort"), length(main.data.names))),
                  latex.addline("Point FE:", rep("yes", length(model.ls))),
                  latex.addline("Country-year FE:", rep("yes", length(model.ls))),
                  latex.addline("Survey FE:", rep(c("yes","yes","--"), each = 1)),
                  latex.addline("Controls:", rep(c("yes","yes","--"), each = 1)),
                  latex.mean.dv(model.ls))
add.lines[[2]][[length(add.lines[[2]])]] <- paste(add.lines[[2]][[length(add.lines[[2]])]], "\\\\ \\hline \\\\[-1.8ex] ")

# ... save main
fileConn <- file(file.path(tab.path, paste0(stub,".tex")))
writeLines(stargazer(model.ls,
                     title="Changes in time to national/regional capital and local development: Alternative weights",
                     keep = treat.vars,
                     multicolumn = F,# se = se,
                     column.labels =collab_w_ruler(column.labels = names(unlist(lapply(main.data.names, get_depvar))), 
                                                   column.separate = c(2,2,1), 
                                                   trim = 5, add.below = NULL),
                     column.separate = c(2,2,1), column.sep.width ="-5pt",
                     dep.var.caption = "",dep.var.labels = rep("", length(model.ls)),
                     covariate.labels = c(names(treat.vars)),
                     font.size = "scriptsize",
                     notes.align = "c", label=stub, align =T,
                     add.lines = add.lines,digits = 3, intercept.top = T,intercept.bottom = F,
                     omit.stat = c("rsq","res.dev","ser"),
                     notes = latex.notes.comb(), notes.label = "", notes.append = F), 
           fileConn)
close(fileConn)

# Save coefs
all.coef.list[[stub]] <- extract_coef(model.ls)



# CONTROL FOR LOCAL ROADS ######
# Not printed in paper, only in Figures 5 and A16
stub <- "roadscontr"

# Prepare
add_road.vars <- function(dn){
  if(dn == "nightlight"){
    list(c(paste0("I(road.dens.",c(1:6)," > 0)")),
         c(paste0("log(road.dens.",c(1:6)," + 1)")))
  } else {
    list(c(paste0("I(inradius.",c(1:6),".km20 > 0)")),
         c(paste0("log(inradius.",c(1:6),".km20 + 1)")))
  }
}

# Estimate
model.ls <- unlist(lapply(main.data.names, function(dn){
  lapply(add_road.vars(dn), function(av){
    form <- make_form(dv = get_depvar(dn), 
                      expl = c(treat.vars, av, 
                               get_contr(dn)), 
                      fe = fe.global, 
                      iv = "0",
                      se = cluster.global )
    felm(form,  data = data.ls[[dn]], keepX = F, keepCX = F)
  })
}), recursive = F)



#... sum of betas
sum.betas <- lapply(model.ls, function(m){
  dm <- deltaMethod(m, paste(treat.vars, collapse = " + "), vcov. = m$clustervcv)
  wt <- lfe::waldtest(m, as.formula(paste("~", paste(treat.vars, collapse = " + "))), type = "cluster")
  dm$pvalue <- wt["p"]
  dm
})

# ... prepare
add.lines <- list(latex.addline_nomc("$\\beta_1 + \\beta_2$:", 
                                     unlist(lapply(sum.betas, function(x){
                                       paste0(round(x[["Estimate"]], 3), 
                                              ifelse(x[["pvalue"]] < .01, "^{***}",ifelse(x[["pvalue"]] < .05, "^{**}",ifelse(x[["pvalue"]] < .1, "^{*}","" ) ) ))
                                     }))),
                  latex.addline_nomc("", unlist(lapply(sum.betas, function(x) paste0("(",round(x[["SE"]], 3), ")")))),
                  latex.addline("Any road by type:", rep(c("yes","no"), length(main.data.names))),
                  latex.addline("Road density by type:", rep(c("no","yes"), length(main.data.names))),
                  latex.addline("Point FE:", rep("yes", length(model.ls))),
                  latex.addline("Country-year FE:", rep("yes", length(model.ls))),
                  latex.addline("Survey FE:", rep(c("yes","yes","--"), each = 2)),
                  latex.addline("Controls:", rep(c("yes","yes","--"), each = 2)),
                  latex.mean.dv(model.ls))
add.lines[[2]][[length(add.lines[[2]])]] <- paste(add.lines[[2]][[length(add.lines[[2]])]], "\\\\ \\hline \\\\[-1.8ex] ")

# ... save main
fileConn <- file(file.path(tab.path, paste0(stub,".tex")))
writeLines(stargazer(model.ls,
                     title="Changes in time to national/regional capital and local development: Controlling for local roads",
                     keep = c(treat.vars),
                     multicolumn = F,# se = se,
                     column.labels = collab_w_ruler(column.labels = names(unlist(lapply(main.data.names, get_depvar))), 
                                                    column.separate = rep(2, length(main.data.names)), 
                                                    trim = 14, add.below = NULL),
                     column.separate = rep(2, length(main.data.names)), 
                     column.sep.width ="-13pt",
                     dep.var.caption = "",dep.var.labels = rep("", length(model.ls)),
                     covariate.labels = c(names(treat.vars)),
                     font.size = "scriptsize",
                     notes.align = "c", label=stub, align =T,
                     add.lines = add.lines,digits = 3, intercept.top = T,intercept.bottom = F,
                     omit.stat = c("rsq","res.dev","ser"),
                     notes = latex.notes.comb(), notes.label = "", notes.append = F), 
           fileConn)
close(fileConn)

# Save coefs
all.coef.list[[stub]] <- extract_coef(model.ls)


# CONTROL FOR POPULATION ##########
# Not printed in paper, only in Figures 5 and A16
stub <- "popcontr"

# ... estimate
model.ls <- lapply(main.data.names, function(dn){
  form <- make_form(dv = get_depvar(dn), expl = c(treat.vars, "log(1+population)", get_contr(dn)), 
                    fe = fe.global, 
                    iv = "0",
                    se = cluster.global )
  m <- felm(form,  data = data.ls[[dn]], keepX = F, keepCX = F)
  m
})


#... sum of betas
sum.betas <- lapply(model.ls, function(m){
  dm <- deltaMethod(m, paste(treat.vars, collapse = " + "), vcov. = m$clustervcv)
  wt <- lfe::waldtest(m, as.formula(paste("~", paste(treat.vars, collapse = " + "))), type = "cluster")
  dm$pvalue <- wt["p"]
  dm
})

# ... prepare
add.lines <- list(latex.addline_nomc("$\\beta_1 + \\beta_2$:", 
                                     unlist(lapply(sum.betas, function(x){
                                       paste0(round(x[["Estimate"]], 3), 
                                              ifelse(x[["pvalue"]] < .01, "^{***}",ifelse(x[["pvalue"]] < .05, "^{**}",ifelse(x[["pvalue"]] < .1, "^{*}","" ) ) ))
                                     }))),
                  latex.addline_nomc("", unlist(lapply(sum.betas, function(x) paste0("(",round(x[["SE"]], 3), ")")))),
                  latex.addline("Point FE:", rep("yes", length(model.ls))),
                  latex.addline("Country-year FE:", rep("yes", length(model.ls))),
                  latex.addline("Survey FE:", rep(c("yes","yes","--"), 1)),
                  latex.addline("Controls:", rep(c("yes","yes","--"), 1)),
                  latex.mean.dv(model.ls))
add.lines[[2]][[length(add.lines[[2]])]] <- paste(add.lines[[2]][[length(add.lines[[2]])]], "\\\\ \\hline \\\\[-1.8ex] ")

# ... save main
fileConn <- file(file.path(tab.path, paste0(stub,".tex")))
writeLines(stargazer(model.ls,
                     title="Changes in time to national/regional capital and local development: Controlling for population",
                     keep = c(treat.vars, "population"),
                     multicolumn = F,# se = se,
                     column.labels = names(unlist(lapply(main.data.names, get_depvar))),
                     #column.separate = c(3,3), column.sep.width ="-10pt",
                     dep.var.caption = "",dep.var.labels = rep("", length(model.ls)),
                     covariate.labels = c(names(treat.vars), "Population (log)"),
                     font.size = "scriptsize",
                     notes.align = "c", label=stub, align =T,
                     add.lines = add.lines,digits = 3, intercept.top = T,intercept.bottom = F,
                     omit.stat = c("rsq","res.dev","ser"),
                     notes = latex.notes.comb(), notes.label = "", notes.append = F), 
           fileConn)
close(fileConn)

# Save coefs
all.coef.list[[stub]] <- extract_coef(model.ls)


# 1966 ROAD NETWORK ##############
# Not printed in paper, only in Figures 5 and A16
stub <- "1966roadnw"


# Estimate
model.ls <- lapply(main.data.names, function(dn){
  y = 1966
  form <- make_form(dv = get_depvar(dn), expl = c(paste0("log(1+",c("natcap.", "regcap."), y, ")"),
                                                 get_contr(dn)), 
                    fe = fe.global, 
                    iv = "0",
                    se = cluster.global )
  felm(form,  data =  data.ls[[dn]], keepX = F, keepCX = F)
})



# ... prepare
add.lines <- list(latex.addline("Point FE:", rep("yes", length(model.ls))),
                  latex.addline("Country-year FE:", rep("yes", length(model.ls))),
                  latex.addline("Survey FE:", rep(c("yes","yes","--"), 1)),
                  latex.addline("Controls:", rep(c("yes","yes","--"), 1)),
                  latex.mean.dv(model.ls))

# ... save main
fileConn <- file(file.path(tab.path, paste0(stub,".tex")))
writeLines(stargazer(model.ls,
                     title="Changes in time to national/regional capital and local development: Stable 1966 road network",
                     keep = treat.vars,
                     multicolumn = F,# se = se,
                     column.labels = names(unlist(lapply(main.data.names, get_depvar))),
                     #column.separate = c(3,3), column.sep.width ="-10pt",
                     dep.var.caption = "",dep.var.labels = rep("", length(model.ls)),
                     covariate.labels = names(treat.vars),
                     font.size = "scriptsize",
                     notes.align = "c", label=stub, align =T,
                     add.lines = add.lines,digits = 3, intercept.top = T,intercept.bottom = F,
                     omit.stat = c("rsq","res.dev","ser"),
                     notes = latex.notes.comb(.85), notes.label = "", notes.append = F), 
           fileConn)
close(fileConn)

# Save coefs
all.coef.list[[stub]] <- extract_coef(model.ls)

##################################
# MAKE ROBUSTNESS CHECK PLOTS ####
# FIGURES 5 and A16 

## Save 
saveRDS(all.coef.list, "data/tempresults.rds")

## Make Plot
source("code/robcheck_plot.R")



# FIGURE A16: BY ROAD YEAR ###################
stub <- "byroadyear"


if(ncore > 1){
  ##### PARALLELIZATION ########
  # Setup cluster
  # ... subset data to what is needed
  dec.data.ls <- lapply(main.data.names, function(dn){
    form <- make_form(dv = get_depvar(dn), expl = c(treat.vars, 
                                                    get_contr(dn)), 
                      fe = fe.global, 
                      iv = "0",
                      se = cluster.global )
    this.data <- data.ls[[dn]]
    na.omit(this.data[,unique(c("year", "cowcode", all.vars(form), 
                                colnames(this.data)[grepl(paste(1945:2017, collapse = "|"), 
                                                          colnames(this.data))]))])
  })
  names(dec.data.ls) <- main.data.names
  
  # ... prepare cluster
  print(paste(Sys.time(), "Setup cluster, load objects"))
  try(registerDoSEQ())
  cl <- makeCluster(getOption("cl.cores", ncore))
  clusterExport(cl, c("get_contr", "get_depvar", "make_form", 
                      "fe.global", "cluster.global", "treat.vars"),
                envir = environment())
  registerDoParallel(cl)
  
  
  # Individual years with road network data
  road.years <- na.omit(as.numeric(unique(gsub("[^\\d]+", "", 
                                               colnames(data.ls[[1]])[grep("regcap.1|regcap.2", 
                                                                           colnames(data.ls[[1]]))], perl = T))))
  road.years <- road.years[order(road.years)]
  
  
  # Estimate
  model.ls <- lapply(main.data.names, function(dn){
    print(dn)
    # Export data to cluster
    this.data <- dec.data.ls[[dn]]
    export.ls <- c("this.data")
    clusterExport(cl, export.ls,
                  envir = environment())
    
    # Estimate
    foreach(y = road.years, .noexport = export.ls, .packages = c("lfe")) %dopar% {
      form <- make_form(dv = get_depvar(dn), expl = c(paste0("log(1+",c("natcap.", "regcap."), y, ")"),
                                                      get_contr(dn)), 
                        fe = fe.global, 
                        iv = "0",
                        se = cluster.global )
      if(all(paste0(c("natcap.", "regcap."), y) %in% colnames(this.data))){
        felm(form,  data = this.data, keepX = F, keepCX = F)
      } else {
        NULL
      }
    }
  })
} else {
  ########## SEQUENTIAL ESTIMATION ##############
  
  # Individual years with road network data
  road.years <- na.omit(as.numeric(unique(gsub("[^\\d]+", "", 
                                               colnames(data.ls[[1]])[grep("regcap.1|regcap.2", 
                                                                           colnames(data.ls[[1]]))], perl = T))))
  road.years <- road.years[order(road.years)]
  
  
  # Estimate
  model.ls <- lapply(main.data.names, function(dn){
    # Estimate
    lapply(road.years, function(y){
      form <- make_form(dv = get_depvar(dn), expl = c(paste0("log(1+",c("natcap.", "regcap."), y, ")"),
                                                      get_contr(dn)), 
                        fe = fe.global, 
                        iv = "0",
                        se = cluster.global )
      if(all(paste0(c("natcap.", "regcap."), y) %in% colnames(data.ls[[dn]]))){
        felm(form,  data = data.ls[[dn]], keepX = F, keepCX = F)
      } else {
        NULL
      }
    })
  })
}

  

# Get Coefs
coef.df <- do.call(rbind, lapply(1:length(main.data.names), function(dn){
  do.call(rbind, lapply(1:length(road.years), function(i){
    m <- model.ls[[dn]][[i]]
    if(class(m) == "try-error" | is.null(m)){
      NULL
    } else {
      pos <- which(grepl("natcap|regcap", row.names(m$coefficients)))
      data.frame(outcome = rep(main.data.names[dn], 2),
                 roadyear = rep(road.years[i], 2),
                 var = treat.vars,
                 coef = m$coefficients[pos, 1],
                 se = diag(m$clustervcv[pos,pos])^.5,
                 stringsAsFactors = F)
    }
  }))
}))

# ... plot
plot.df <- coef.df
plot.df$outcome.name <- factor(names(unlist(lapply(plot.df$outcome, get_depvar))), 
                               levels = names(unlist(lapply(main.data.names, get_depvar))))
plot.df$var.name <- ifelse(plot.df$var == treat.vars[1], names(treat.vars)[1], names(treat.vars)[2])


# ... print
png(file.path(fig.path, paste0("figa16_byroadyr.png")), width = 8, height = 4, units = "in", res = 300)
p <- ggplot(plot.df, aes(x = roadyear, y = coef)) +
  geom_hline(yintercept = 0, lty = 2, col = "grey") +
  geom_point() +
  geom_segment(aes(y = coef + se * 1.96,
                   yend = coef - se * 1.96,
                   x = roadyear, xend = roadyear)) +
  facet_wrap(var.name ~ outcome.name , scales = "free_y") +
  xlab("Road network year") + 
  ylab("Marginal effect of time to\nnat./reg. capital") +
  scale_x_continuous(breaks = seq(1970, 2010, by = 20)) +
  theme_minimal()
print(p)
dev.off()


# Stop Cluster
if(ncore > 1){
  stopCluster(cl)
  rm(cl)
}




# FIGURE A17: NIGHLIGHTS: DIFFERENT UNIT SIZES ###
stub <- "nlunitsize"


# Prepare
outcomes <- c("nightlight.any", "nightlight.log", "nightlight.pc.log")
names(outcomes) <- c("Any Light (0/100)", "Light (log)", "Light/capita (log)")

# ... unit sizes
unit.sizes <- c(100, 200, 400, 800, 1600, 3200, 6400)

# Estimate // loads and discards data on the fly to save RAM
model.ls <- lapply(seq_along(unit.sizes), function(i){
  print(i)
  data <- readRDS(file.path("data", 
                            paste0("data_nightlight_",unit.sizes[i], ".rds")))
  form <- make_form(dv = get_depvar("nightlight"), 
                    expl = c(treat.vars,get_contr("nightlight")), 
                    fe = fe.global, 
                    iv = "0",
                    se = cluster.global )
  felm(form,  data = data, keepX = F, keepCX = F)
})


# Get Coefs
coef.df <- do.call(rbind, lapply(seq_along(unit.sizes), function(i){
  m <- model.ls[[i]]
  if(class(m) == "try-error"){
    NULL
  } else {
    data.frame(outcome = rep(names(get_depvar("nightlight")), length(treat.vars)),
               size = rep(i, length(treat.vars)),
               var = treat.vars,
               coef = m$coefficients[treat.vars, 1]/abs(mean(m$response)),
               se = (diag(m$clustervcv[treat.vars,treat.vars])^.5)/abs(mean(m$response)),
               stringsAsFactors = F)
  }
}))

# ... plot
plot.df <- coef.df
plot.df$var.name <- ifelse(plot.df$var == treat.vars[1], 
                           names(treat.vars)[1], names(treat.vars)[2])


# ... make figure
png(file.path(fig.path, paste0("figa17_nightlight_unitsize_main.png")), 
    width = 6, height = 4, units = "in", res = 300)
p <- ggplot(plot.df, aes(x = size, y = coef)) +
  geom_hline(yintercept = 0, lty = 2, col = "grey") +
  geom_point() +
  geom_segment(aes(y = coef + se * 1.96,
                   yend = coef - se * 1.96,
                   x = size, xend = size)) +
  facet_wrap(~ var.name , nrow = 1, scales = "free_y") +
  xlab("Unit size in square km") + 
  ylab("Std. marginal effect of time to nat./reg. capital \n on nightlights/capita (log)") +
  scale_x_continuous(breaks = seq_along(unit.sizes),
                     labels = unit.sizes) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5)) +
  theme_minimal()
print(p)
dev.off()

# CLEANUP ######################

# Remove main data
rm(data.ls)





